# Alexander F. Gazmararian
# afg2@princeton.edu

library(here)
library(tidyverse)
library(tidylog)
library(modelsummary)

#load fair data
fair <- readRDS(here("data", "fair_conjoint.rds"))
fair <- subset(fair, !is.na(selected))

#load national data
nat <- readRDS(here("data", "national_conjoint.rds"))

#clean up variable names
fair <-
  fair %>%
  rename(
    Free.Retraining.Program = `Free Retraining Program`,
    Retrained.Worker.Salary = `Retrained Worker Salary`,
    Worker.Retraining.Time = `Worker Retraining Time`,
    Income.Support.During.Retraining = `Income Support During Retraining`,
    Relocation.Support = `Relocation Support`,
    Community.Investment = `Community Investment`,
    Benefit.Support = `Benefit Support`
  )
fair$weights <- fair$weights2

#specify equation
f <- selected ~ Community.Investment + Free.Retraining.Program + Income.Support.During.Retraining + 
  Worker.Retraining.Time + Benefit.Support + Retrained.Worker.Salary + Relocation.Support

#estimate models
m.nat <- lm(f, nat, weights = weights)
summary(m.nat)
m.fair <- lm(f, fair, weights = weights)
summary(m.fair)

#extract and combine coefficients
cj <- modelplot(list(m.nat, m.fair), vcov = ~ ResponseId, draw = FALSE)
#combine and estimate difference
fair$sample <- "Fossil Fuel"
nat$sample <- "National"
g <- rbind(
  subset(fair,select=c(sample,selected,ResponseId,task,Free.Retraining.Program:Community.Investment,weights)),
         subset(nat,select=c(sample,selected,ResponseId,task,Free.Retraining.Program:Community.Investment,weights))
  )
glimpse(g)
m.int <- lm(update(f, ~ . * sample), g, weights = weights)
#Create table
tab <- here("output", "tables", "si_tab_8.2.txt")
modelsummary(
  m.int,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = c(
    "(Intercept)"="Intercept",
    
    "Community.InvestmentBroadband"="Community Investment: Broadband",
    "Community.InvestmentBroadband:sampleNational" = "Community Investment: Broadband $\\times$ National",
    
    "Community.InvestmentHousing for new residents" = "Community Investment: Housing for new residents",
    "Community.InvestmentHousing for new residents:sampleNational" = "Community Investment: Housing for new residents $\\times$ National",
    
    "Free.Retraining.ProgramFor jobs in clean energy" = "Free Retraining Program: Clean Energy",
    "Free.Retraining.ProgramFor jobs in clean energy:sampleNational" = "Free Retraining Program: Clean Energy $\\times$ National",
    
    "Free.Retraining.ProgramFor jobs in health care"= "Free Retraining Program: Health Care",
    "Free.Retraining.ProgramFor jobs in health care:sampleNational"= "Free Retraining Program: Health Care $\\times$ National",
    
    "Income.Support.During.Retraining$400 per week for workers"="Income Support During Retraining: \\$400 per week for workers",
    "Income.Support.During.Retraining$400 per week for workers:sampleNational"="Income Support During Retraining: \\$400 per week for workers $\\times$ National",
    
    "Income.Support.During.RetrainingNone"="Income Support During Retraining",
    "Income.Support.During.RetrainingNone:sampleNational"="Income Support During Retraining $\\times$ National",
    
    "Worker.Retraining.Time2 years"="Worker Retraining Time: 2 years",
    "Worker.Retraining.Time2 years:sampleNational"="Worker Retraining Time: 2 years $\\times$ National",
    
    "Worker.Retraining.Time3-6 months"="Worker Retraining Time: 3-6 Months",
    "Worker.Retraining.Time3-6 months:sampleNational"="Worker Retraining Time: 3-6 Months $\\times$ National",
    
    "Benefit.SupportFunding for worker health care"="Benefit Support: Worker Health Care",
    "Benefit.SupportFunding for worker health care:sampleNational"="Benefit Support: Worker Health Care $\\times$ National",
    
    "Benefit.SupportFunding for worker pensions"="Benefit Support: Pensions",
    "Benefit.SupportFunding for worker pensions:sampleNational"="Benefit Support: Pensions $\\times$ National",
    
    "Retrained.Worker.Salary$100,000 per year"="Retrained Worker Salary: \\$100,000 per year",
    "Retrained.Worker.Salary$100,000 per year:sampleNational"="Retrained Worker Salary: \\$100,000 per year $\\times$ National",
    
    "Retrained.Worker.Salary$50,000 per year"="Retrained Worker Salary: \\$50,000 per year",
    "Retrained.Worker.Salary$50,000 per year:sampleNational"="Retrained Worker Salary: \\$50,000 per year $\\times$ National",
    
    "Relocation.SupportVoucher to cover moving expenses"="Relocation Support: Voucher to Cover Moving Expenses",
    "Relocation.SupportVoucher to cover moving expenses:sampleNational"="Relocation Support: Voucher to Cover Moving Expenses $\\times$ National"
  ),
  vcov = ~ ResponseId,
  gof_omit = "IC|Log|RMSE|Std",
  escape = FALSE,
  output = "latex"
) %>%
  cat(., file = tab)
tab2 <- readLines(tab)
cat(tab2[-c(1,2,66,68)], file = tab)



#Create figure
df.int <- modelplot(m.int, vcov = ~ ResponseId, draw = FALSE)
df.int %>%
  filter(grepl(":sampleNational",term)) %>%
  ggplot() +
  geom_hline(yintercept=0, color = "grey") +
  geom_pointrange(aes(x=term,y=estimate,ymin=conf.low,ymax=conf.high),
                  position = position_dodge(.5)) +
  theme_bw() +
  coord_flip() +
  theme(
    panel.grid = element_blank()
  )
#add together and get p-value for difference
df.int.merge <- subset(df.int, select = c(p.value,term), grepl(":sample", term))
names(df.int.merge)[1] <- "p.df"
df.int.merge$term <- gsub(":sampleNational", "", df.int.merge$term)
plot.df <-
  cj %>%
  left_join(., df.int.merge, by = "term") %>%
  mutate(p.df = ifelse(model=="(1)", NA, p.df),
         p.df = case_when(
           p.df >= .05 & p.df <.1 ~ "*",
           p.df < 0.05 & p.df >= .01 ~ "**",
           p.df < 0.01 ~ "***",
           T ~ NA
         )) %>%
  filter(term!="(Intercept)" & !grepl("task|ResponseId",term)) %>%
  mutate(model=ifelse(model=="(1)","National", "Fossil Fuel"),
         lb90 = estimate + std.error * qnorm(.05),
         ub90 = estimate + std.error * qnorm(.95))
plot.df <-
  plot.df %>%
  mutate(
    attribute = case_when(
      grepl("Benefit.Support", term) ~ "Benefit Support",
      grepl("Community.Investment", term) ~ "Community Investment",
      grepl("Free.Retraining", term) ~ "Free Retraining Program",
      grepl("Relocation", term) ~ "Relocation Support",
      grepl("Retrained.Worker", term) ~ "Retrained Worker Salary",
      grepl("Worker.Retraining.Time", term) ~ "Worker Retraining Time",
      grepl("Income.Support", term) ~ "Income Support During Retraining"
    )
  )
#relevel
baseline <- cbind(estimate = 0, std.error = 0)
reflevels <- c(
  levels(g$Benefit.Support)[[1]],
  levels(g$Community.Investment)[[1]],
  levels(g$Free.Retraining)[[1]],
  levels(g$Relocation)[[1]],
  levels(g$Worker.Retraining.Time)[[1]],
  levels(g$Income.Support)[[1]]
)
attributelevels <- c("Benefit Support", "Community Investment", "Free Retraining Program", "Relocation Support", "Retrained Worker Salary", "Income Support During Retraining")
#baseline dataframe
base.df <- data.frame(term = reflevels, baseline, attribute = attributelevels)
base.df <- mutate(base.df, across(estimate:std.error, ~ as.numeric(.x)))
# Merge
con <- bind_rows(plot.df, base.df)
# Relevel factors
con$term[con$term == "None" & con$Attribute == "Benefit Support"] <- "No Benefit Support"
con$term[con$term == "None" & con$Attribute == "Income Support During Retraining"] <- "No Income Support"
con$term[con$term == "None" & con$Attribute == "Relocation Support"] <- "No Relocation Support"
con$term[con$term == "$200 per week for all community members"] <- "$200 per week for all"
#Rename term
con <-
  con %>%
  mutate(
    term = case_when(
      term == "Community.InvestmentBroadband" ~ "Broadband",
      term == "Community.InvestmentHousing for new residents" ~ "Housing for new residents",
      term == "Free.Retraining.ProgramFor jobs in clean energy" ~ "For jobs in clean energy",
      term == "Free.Retraining.ProgramFor jobs in health care" ~ "For jobs in health care",
      term == "Income.Support.During.Retraining$400 per week for workers" ~ "$400 per week for workers",
      term == "Income.Support.During.RetrainingNone" ~ "None",
      term == "Worker.Retraining.Time2 years" ~ "2 years",
      term == "Worker.Retraining.Time3-6 months" ~ "3-6 months",
      term == "Benefit.SupportFunding for worker health care" ~ "Funding for worker health care",
      term == "Benefit.SupportFunding for worker pensions" ~ "Funding for worker pensions",
      term == "Relocation.SupportVoucher to cover moving expenses" ~ "Voucher to cover moving expenses",
      term == "Retrained.Worker.Salary$100,000 per year" ~ "$100,000 per year",
      term == "Retrained.Worker.Salary$50,000 per year" ~ "$50,000 per year",
      T ~ NA_character_
    )
  )
#Create data frames
benefit <- data.frame(attribute = "Benefit Support", term = "Benefit Support")
community <- data.frame(attribute = "Community Investment", term = "Community Investment")
train <- data.frame(attribute = "Free Retraining Program", term = "Free Retraining Program")
income <- data.frame(attribute = "Income Support During Retraining", term = "Income Support During Retraining")
locate <- data.frame(attribute = "Relocation Support", term = "Relocation Support")
salary <- data.frame(attribute = "Retrained Worker Salary", term = "Retrained Worker Salary")
time <- data.frame(attribute = "Worker Retraining Time", term = "Worker Retraining Time")
con <- bind_rows(con, benefit, community, train, income, locate, salary, time)
# Arrange by level
con <-
  con %>%
  arrange(attribute, !is.na(estimate))
# Fix names
con$term[[1]] <- "**Benefit Support**<br>Baseline: None"
con$term[[7]] <- "**Community Investment**<br>Baseline: Schools"
con$term[[13]] <- "**Free Retraining Program**<br>Baseline: For jobs in manufacturing"
con$term[[19]] <- "**Income Support**<br>Baseline: $200 per week for all"
con$term[[25]] <- "**Relocation Support**<br>Baseline: None"
con$term[[29]] <- "**Retrained Worker Salary**<br>Baseline: $75,000 per year"
con$term[[35]] <- "**Worker Retraining Time**<br>Baseline: 1 year"

con <- subset(con, !is.na(term))
glimpse(con)

# Create levels
con$term <- factor(con$term, ordered = TRUE, levels = rev(c(
  "**Relocation Support**<br>Baseline: None", "Voucher to cover moving expenses",
  "**Community Investment**<br>Baseline: Schools", "Broadband", "Housing for new residents",
   "**Benefit Support**<br>Baseline: None", "Funding for worker health care", "Funding for worker pensions",
   "**Income Support**<br>Baseline: $200 per week for all", "$400 per week for workers",  "None",
  "**Free Retraining Program**<br>Baseline: For jobs in manufacturing", "For jobs in clean energy", "For jobs in health care",
  "**Retrained Worker Salary**<br>Baseline: $75,000 per year", "$100,000 per year",  "$50,000 per year",
  "**Worker Retraining Time**<br>Baseline: 1 year", "2 years", "3-6 months"
)))

back_line <- 0.25
pos_dodge <- .7
plot.cj <-
  con %>%
  filter(!is.na(term)) %>%
  ggplot(aes(x=term,y=estimate,color=model,shape=model)) +
  geom_vline(xintercept = "**Community Investment**<br>Baseline: Schools", color = "gray", size = back_line) +
  geom_vline(xintercept = "**Benefit Support**<br>Baseline: None", color = "gray", size = back_line) +
  geom_vline(xintercept = "**Free Retraining Program**<br>Baseline: For jobs in manufacturing", color = "gray", size = back_line) +
  geom_vline(xintercept = "**Income Support**<br>Baseline: $200 per week for all", color = "gray", size = back_line) +
  geom_vline(xintercept = "**Retrained Worker Salary**<br>Baseline: $75,000 per year", color = "gray", size = back_line) +
  geom_vline(xintercept = "**Worker Retraining Time**<br>Baseline: 1 year", color = "gray", size = back_line) +
  geom_hline(yintercept=0, color = "grey") +
  geom_errorbar(aes(ymin=conf.low,ymax=conf.high),width=0,position=position_dodge(pos_dodge)) +
  geom_errorbar(aes(ymin=lb90,ymax=ub90),width=0,size=1.25,position=position_dodge(pos_dodge)) +
  geom_point(size = 3,fill="white",position=position_dodge(pos_dodge)) +
  geom_text(aes(x=term,label=p.df), size = 6, show.legend = FALSE) +
  scale_color_grey(na.translate = FALSE, start = 0, end = .5) +
  scale_shape_manual(values = c(21,24), na.translate = FALSE) +
  theme_bw(base_size = 14) +
  coord_flip() +
  scale_y_continuous(labels = scales::percent) +
  labs(color="",shape="",y="Change in Probability of Support", x="") +
  theme(
    panel.grid = element_blank(),
    axis.text.y = ggtext::element_markdown(),
    legend.position = c(.85,.08),
    legend.background = element_blank()
  )
plot.cj
ggsave(
  plot.cj,
  filename = here("output", "figures", "fig_1.png"),
  scale = 1.25,
  width = 6.5,
  height = 6.25,
  dpi = 300
)
